Spectra Story

# Libraries
library(tidyverse)
library(gmRi)
library(rnaturalearth)
library(scales)
library(sf)
library(gt)
library(patchwork)
library(lmerTest)
library(emmeans)
library(merTools)
library(tidyquant)
library(ggeffects)
library(performance)
library(gtsummary)
library(gt)
library(sizeSpectra)
library(ggdist)
library(pander)
library(ggh4x)
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")
conflicted::conflicts_prefer(lmerTest::lmer)

# Processing functions
source(here::here("Code/R/Processing_Functions.R"))

# ggplot theme
theme_set(
  theme_gmri(
    rect = element_rect(fill = "white", color = NA), 
    strip.text.y = element_text(angle  = 0),
    axis.text.x = element_text(size = 8),
    legend.position = "bottom"))



# vectors for factor levels
area_levels <- c("GoM", "GB", "SNE", "MAB")
area_levels_long <- c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")
area_levels_all <- c("Northeast Shelf", area_levels)
area_levels_long_all <- c("Northeast Shelf", area_levels_long)

# table to join for swapping shorthand for long-hand names
area_df <- data.frame(
  area = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"),
  survey_area = c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),
  area_titles = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))


# Degree symbol
deg_c <- "\u00b0C"



# Function to get the Predictions
# Drop effect fits that are non-significant  ###
get_preds_and_trendsignif <- function(mod_x){
  modx_preds <- as.data.frame(
    # Model predictions
    ggpredict(
      mod_x, 
      terms = ~ est_year * survey_area * season) ) %>% 
    mutate(
      survey_area = factor(group, levels = area_levels),
      season = factor(facet, levels = c("Spring", "Fall")))
  
    # Just survey area and year
    modx_emtrend <- emtrends(
      object = mod_x,
      specs =  ~ survey_area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
  
    # Preds with signif
    modx_preds %>% left_join(select(modx_emtrend, survey_area, season, non_zero))
  
}

# Broad Distribution
#log2 bins for easy-of-access
#' @title Build Log 2 Bin Structure Dataframe
#' 
#' @description Used to build a dataframe containing equally spaced log2 bins for
#' size spectra analysis. Contains details on the left and right limits, midpoint, bin width, 
#' and a text label for the bins. log2bin number ascends with increasing size for eeasy plotting.
#'
#' @param log2_min 
#' @param log2_max 
#' @param log2_increment 
#'
#' @return
#' @export
#'
#' @examples
define_log2_bins <- function(log2_min = 0, log2_max = 13, log2_increment = 1){
  
  # How many bins
  n_bins  <- length(seq(log2_max, log2_min + log2_increment, by = -log2_increment))
  
  # Build Equally spaced log2 bin df
  log2_bin_structure <- data.frame(
    "log2_bins" = as.character(seq(n_bins, 1, by = -1)),
    "left_lim"  = seq(log2_max - log2_increment, log2_min, by = -log2_increment),
    "right_lim" = seq(log2_max, log2_min + log2_increment, by = -log2_increment)) %>% 
    mutate(
      bin_label    = str_c(round(2^left_lim, 3), " - ", round(2^right_lim, 3), "g"),
      bin_width    = 2^right_lim - 2^left_lim,
      bin_midpoint = (2^right_lim + 2^left_lim) / 2) %>% 
    arrange(left_lim)
  
  return(log2_bin_structure)
}






#' @title Assign Manual log2 Bodymass Bins - By weight
#'
#' @description Manually assign log2 bins based on individual length-weight bodymass 
#' in increments of 1 on the log2 scale. Returns data with bins assigned based on individual
#' length-weight biomass
#' 
#' Uses maximum weight, and its corresponding bin as the limit.
#'
#' @param wmin_grams Catch data prepared for mle calculation, use prep_wmin_wmax
#' @param log2_increment Equally spaced increments to use for log 2 bin sizes. Default = 0.5.
#'
#' @return
#' @export
#'
#' @examples
assign_log2_bins <- function(wmin_grams, log2_increment = 1){
  
  
  #### 1. Set up bodymass bins
  
  # filter missing weights
  size_data <- wmin_grams %>% 
    filter(wmin_g > 0,
           is.na(wmin_g) == FALSE,
           wmax_g > 0,
           is.na(wmax_g) == FALSE)
  
  # Get bodymass on log2() scale
  size_data$log2_weight <- log2(size_data$ind_weight_g)
  
  # Set up the bins - Pull min and max weights from data available
  #min_bin <- floor(min(size_data$log2_weight))
  min_bin <- 0
  max_bin <- ceiling(max(size_data$log2_weight))
  
  
  # Build a bin key, could be used to clean up the incremental assignment or for apply style functions
  log2_bin_structure <- define_log2_bins(
    log2_min = min_bin, 
    log2_max = max_bin, 
    log2_increment = log2_increment)
  
  
  
  # Loop through bins to assign the bin details to the data
  log2_assigned <- log2_bin_structure %>%
    split(.$log2_bins) %>%
    map_dfr(function(log2_bin){
      
      # limits and labels
      l_lim   <- log2_bin$left_lim
      r_lim   <- log2_bin$right_lim
      bin_num <- as.character(log2_bin$log2_bin)
      
      # assign the label to the appropriate bodymasses
      size_data %>% mutate(
        log2_bins = ifelse( between(log2_weight, l_lim, r_lim), bin_num, NA),
        log2_bins = as.character(log2_bins)) %>%
        drop_na(log2_bins)
      
    })
  
  # Join in the size bins
  log2_assigned <- left_join(log2_assigned, log2_bin_structure, by = "log2_bins")
  
  # return the data with the bins assigned
  return(log2_assigned)
  
}

Storycrafting - Northeast Shelf Spectra

I think we need to step back a bit and find the main story again–which is about how size structure of the NES fish community has changed over time, whether that varies by region, and what covariates are associated with those changes.
The Aug draft you developed really focused on comparing regression models to consider drivers based on the fish community size spectrum as viewed from the perspective of the Wigley species set vs. the full species set.
However, I didn’t see figures of those spectra in that draft (maybe there are somewhere else and I’ve lost track). The recent slide deck really focuses on seasonal differences in the size spectra and influence of temperature. - KMills

This document will be split into two sections (the second potentially moving to its own doc).

Part 1: Summary of Trends

For the first section the focus is on understanding changes in the size structure. This will be done using three metrics: 1. length and bodymass spectra 2. large and small fish indices 3. Body length/mass size quantiles changes over time.

This section should address: Has the size distribution changed over time, where has it changed, and some nuance around what aspects about the distribution have shifted.

Starting from Square Uno

I’d like for us to start by rebuilding the story of how the size structure has changed over time and regional differences.
As such, I would like for us to look at plots of the size spectrum slope for (1) the all-species length spectrum and (2) the Wigley species biomass spectrum from annual, seasonal and regional perspectives.
I am thinking this would be two panel plots (all-spp length, Wigley-spp biomass) with 5 subpanels each (NES, GOM, GB, SNE, MAB) that are set up to show lines by seasonal and annual time steps, much like the plot on the right of slide 7 in your most recent slide deck…except with NES added and an annual line added. Keeping the significant trend lines in is helpful. This is the first step, so when you have that ready, please share it with the group. - KMills

# Data used for Wigley estimates
wigley_in <- read_csv(here::here("Data/processed/wigley_species_trawl_data.csv"))
finfish_in <- read_csv(here::here("Data/processed/finfish_trawl_data.csv"))

# Groundfish targeted stocks
targeted_species <- c(
  "atlantic cod",
  "haddock",
  "american plaice",
  "white hake",
  "witch flounder",
  "winter flounder",
  "yellowtail flounder",
  "pollock",
  "acadian redfish",
  "atlantic halibut",
  "windowpane flounder",
  "ocean pout",
  "atlantic wolffish")

# filter the groundfish complex
gfish_comp <- filter(wigley_in, comname %in% targeted_species)

Aside: Biomass Fraction in Each Community

If we need to make a decision about which community we want to use, here is the different biomass totals for each:

# glimpse(wigley_in)

# Load raw
# General tidying only, removal of strata outside our study area, Spring and Fall only
# (removes inshore and rarely/inconsistently sampled strata etc.)
# shrimps removed
trawl_basic <- read_csv(here::here("Data/processed/trawl_clean_all_data.csv"))

# Biomass totals
raw_totals <- trawl_basic %>% 
  distinct(est_towdate, cruise6, station, stratum, tow, comname,  biomass_kg) %>%   pull(biomass_kg) %>% 
  sum()

wig_biomass_total <- wigley_in %>% 
  distinct(est_towdate, cruise6, station, stratum, tow, comname,  biomass_kg) %>%   pull(biomass_kg) %>% 
  sum()

ffish_biomass_total <- finfish_in %>% 
  distinct(est_towdate, cruise6, station, stratum, tow, comname,  biomass_kg) %>%   
  pull(biomass_kg) %>% 
  sum()

gfish_total  <- gfish_comp %>% 
  distinct(est_towdate, cruise6, station, stratum, tow, comname,  biomass_kg) %>%   
  pull(biomass_kg) %>% 
  sum()

# Make a table for display
tibble(
  "Community" = c("All (including crustaceans)", "All finfish", "Wigley Species Subset", "Groundfish Complex"),
  "Biomass Total" = c(raw_totals, ffish_biomass_total, wig_biomass_total, gfish_total)) %>% 
  mutate(`% of Total` = (`Biomass Total` / raw_totals)*100) %>% 
  gt()

Community	Biomass Total	% of Total
All (including crustaceans)	3781548	100.00000
All finfish	3708046	98.05629
Wigley Species Subset	3539719	93.60501
Groundfish Complex	932191	24.65104

The general tradeoffs are this:

The more selective we are about which species we include, then the more targeted we can be about questions of “why” and the more assured about their consistent sampling/selectivity for the gear.

If we are too selective then we may miss broader patterns. If we aren’t selective enough, then we may likely be incorporating unnecessary noise that may never be resolved.

Using the groundfish complex seems safe in the sense that those species are all well sampled and are well studied. They also are most likely to show a predictable response to fisheries landings.

Using the whole finfish community has me feeling enigmatic whether I’m looking at a community resilient to disturbance, or just noise stemming from inclusion of poorly/inconsistently sampled species.

The Groundfish Complex Species

I had the thought that it may be easier to identify the sensitivity to fisheries impacts by estimating spectra for a subset of species that make up the groundfish complex.

These are the species in question:

The 20 large-mesh stocks are:

Targeted Stocks:

1. Gulf of Maine cod
   
2. Georges Bank cod
   
3. Gulf of Maine haddock
   
4. Georges Bank haddock
   
5. Georges Bank yellowtail flounder
   
6. Southern New England/Mid-Atlantic yellowtail flounder
   
7. Cape Cod/Gulf of Maine yellowtail flounder
   
8. pollock
   
9. American plaice
   
10. witch flounder
    
11. white hake
    
12. Georges Bank winter flounder
    
13. Gulf of Maine winter flounder
    
14. Southern New England/Mid-Atlantic winter flounder
    
15. redfish
    
16. Atlantic halibut

Non-Targeted Stocks:

17. Gulf of Maine/Georges Bank (northern) windowpane flounder
    
18. Southern New England/Mid-Atlantic (southern) windowpane flounder
    
19. ocean pout
    
20. Atlantic wolffish

# # # Run MAB with and without spiny dogfish
# # # Run the year, season, area fits for the filtered data
# gf_bmspectra <- group_binspecies_spectra(
#   ss_input = gfish_comp,
#   grouping_vars = c("est_year", "season", "survey_area"),
#   abundance_vals = "numlen_adj",
#   length_vals = "length_cm",
#   use_weight = TRUE,
#   isd_xmin = 1,
#   global_min = TRUE,
#   isd_xmax = NULL,
#   global_max = FALSE,
#   bin_width = 1,
#   vdiff = 2)

# # Save those
# write_csv(
#   gf_bmspectra,
#   here::here("Data/processed/groundfish_complex_bmspectra.csv"))

A. Size Spectra Trends

For the figures in part 1: non-zero trends over time have been highlighted in the appropriate plot panels with linear trends overlaid.

# Load the three regional datasets
# wigley_lenspectra <- read_csv(
#   here::here("Data/processed/wigley_species_length_spectra.csv"))  %>% 
#   mutate(survey_area = factor(survey_area, levels = area_levels))
wigley_bmspectra <- read_csv(
  here::here("Data/processed/wigley_species_bodymass_spectra.csv")) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels))
finfish_lenspectra <- read_csv(
  here::here("Data/processed/finfish_length_spectra.csv")) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels))

# Load Shelf
finfish_lenspectra_shelf <- read_csv(
  here::here("Data/processed/shelfwide_finfish_length_spectra.csv"))
wigley_bmspectra_shelf <- read_csv(
  here::here("Data/processed/shelfwide_wigley_species_bodymass_spectra.csv"))




# Join them together
wigley_bmspectra_all <- bind_rows(wigley_bmspectra, wigley_bmspectra_shelf) %>% 
  mutate(metric = "bodymass_spectra",
         community = "Wigley Subset",
         survey_area = factor(survey_area, levels = area_levels_all))
ffish_lenspectra_all <- bind_rows(finfish_lenspectra, finfish_lenspectra_shelf) %>% 
    mutate(metric = "length_spectra",
           community = "Finfish Community",
         survey_area = factor(survey_area, levels = area_levels_all))
# both metrics+communities
spectra_results <-  bind_rows(wigley_bmspectra_all, ffish_lenspectra_all) %>% 
  mutate(
    survey_area = fct_relevel(survey_area, area_levels_all),
    metric = fct_rev(metric))

Wigley Species - Bodymass Spectra

# Model linear trends in time
bmspectra_mod_wig_all <- lmer(
  formula = b ~ scale(est_year) * survey_area * season + (1|est_year),
  data = wigley_bmspectra_all)

# Get the predictions and flag whether they are significant
bmspectra_trend_predictions <- get_preds_and_trendsignif(bmspectra_mod_wig_all) %>% 
  mutate(model_id = "bodymass_spectra-Wigley Subset") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric),
         survey_area = factor(survey_area, levels = area_levels_all))





# Plot over observed data
# Contrast seasonal differences




# Make the plot
bmspectra_trend_p <- ggplot() +
   geom_ribbon(
    data = filter(bmspectra_trend_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    data = wigley_bmspectra_all,
    aes(est_year, b, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    data = wigley_bmspectra_all,
    aes(est_year, b, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  geom_line(
    data = filter(bmspectra_trend_predictions, non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(survey_area ~ metric*community, scales = "free") +
  labs(
    x = "Year",
    y = "Exponent of Body Mass Distribution",
    linetype = "Linetype",
    color = "Season",
    fill = "Season")
 
# Show plot
bmspectra_trend_p

# # Summary tables
# bmspectra_mod_wig_all %>% 
#   gtsummary::tbl_regression()   %>% 
#   add_global_p(keep = T) %>% 
#   bold_p(t = 0.05)  %>%
#   bold_labels() %>%
#   italicize_levels() %>% 
#   as_gt() %>% 
#   tab_header(title = "Wigley Community - Mass Spectra Trends")

Full Finfish Community - Length Spectra

# Model linear trends in time
lenspectra_mod_all <- lmer(
  formula = b ~ scale(est_year) * survey_area * season + (1|est_year),
  data = ffish_lenspectra_all)

# Get the predictions and flag whether they are significant
lenspectra_trend_predictions <- get_preds_and_trendsignif(lenspectra_mod_all) %>% 
  mutate(model_id = "length_spectra-Finfish Community") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric),
         survey_area = factor(survey_area, levels = area_levels_all))





# Plot over observed data
# Contrast seasonal differences




# Make the plot
lenspectra_trend_p <- ggplot() +
   geom_ribbon(
    data = filter(lenspectra_trend_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    data = ffish_lenspectra_all,
    aes(est_year, b, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    data = ffish_lenspectra_all,
    aes(est_year, b, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  geom_line(
    data = filter(lenspectra_trend_predictions, non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(survey_area ~ metric*community, scales = "free") +
  labs(
    x = "Year",
    y = "Exponent of Body Length Distribution",
    linetype = "Linetype",
    color = "Season",
    fill = "Season")
 
# Show the plot
lenspectra_trend_p

# # Summary tables
# lenspectra_mod_all  %>% 
#   gtsummary::tbl_regression()   %>% 
#   add_global_p(keep = T) %>% 
#   bold_p(t = 0.05)  %>%
#   bold_labels() %>%
#   italicize_levels() %>% 
#   as_gt() %>% 
#   tab_header(title = "Full Finfish Community - Length Spectra Trends")

Both Community Size Spectra

# Left side:
# All species  species length spectrum 
# color = season vs. annual
# facet_wrap(~survey_area)

# Right side:
# Wigley species biomass spectrum
# color = season vs. annual
# facet_wrap(~survey_area)


# Plot over observed data
# Contrast seasonal differences
(
  (lenspectra_trend_p + theme(strip.text.y = element_blank()) | 
  bmspectra_trend_p)
  ) + plot_layout(guides = "collect")

Groundfish Complex Spectra

A couple surprises here.

I was anticipating that these would show a more clear period of change in Northern regions when fisheries for these species were crashing. I am surprised that this is not really the case.

I also anticipated these to show a less powerful seasonal difference. These species are more bottom-associated and likely to have more consistent catchability from a trawl survey. Seeing a much tighter grouping between the seasons AND seeing a semblance of temporal autocorrelation give me a sense that we are following a community through time, and not catching random noise based on what is in the area that time of year and may not be there next year.

(details on the species below)

# Read them in, do some plotting
gf_bmspectra <- read_csv(
  here::here("Data/processed/groundfish_complex_bmspectra.csv")) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels)) %>% 
  mutate(est_year = as.numeric(est_year))


# Join them together
gf_results <- bind_rows(
   list(
     "Wigley Full" = wigley_bmspectra,
     "Groundfish Complex" =  gf_bmspectra), 
   .id = "community") %>% 
  mutate(metric = "bodymass_spectra",
         survey_area = fct_relevel(survey_area, area_levels))


# Get trends for groundfish complex
gfish_mod <- lmer(
  formula = b ~ scale(est_year) * survey_area * season + (1|est_year),
  data = gf_bmspectra)
# check_model(gfish_mod)

#  Model for wigley species
bmspectra_mod_wig <- lmer(
  formula = b ~ scale(est_year) * survey_area * season + (1|est_year),
  data = wigley_bmspectra)

# Put predictions
gf_sens_predictions <- bind_rows(list(
  "bodymass_spectra-Wigley Full" = get_preds_and_trendsignif(bmspectra_mod_wig),
  "bodymass_spectra-Groundfish Complex" = get_preds_and_trendsignif(gfish_mod)), 
  .id = "model_id") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric))



# Plot the differences
ggplot() +
  geom_ribbon(
    data = filter(gf_sens_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    data = gf_results,
    aes(est_year, b, color = season),
    size = 0.8) +
  geom_line(
    data = filter(gf_sens_predictions, non_zero == TRUE),
    aes(x, predicted, color = season),
    linewidth = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(survey_area ~ metric*community, scales = "free") +
  labs(
    x = "Year",
    y = "Exponent of Size Spectra",
    title = "Comparing Wigley Community to Groundfish Complex",
    subtitle = "Spring = Dogfish in MAB, Fall = Dogfish in GB")

B. Median Length & Weight Trends

I also think it would be valuable to build out the story of size change further, so am thinking that plots like those of average (or median) length (all species) and weight (wigley species) for NES and the sub-regions, like those on p. 12 of the attached file above (the very first draft from late 2022) would be useful. - KMills

# Load the median weight/length data
wigley_size_df <- read_csv(here::here("Data/processed/wigley_species_size_summary.csv"))
finfish_size_df <- read_csv(here::here("Data/processed/finfish_species_length_summary.csv"))

# shelf-wide summary
finfish_lengths_shelf <- read_csv(here::here("Data/processed/shelfwide_finfish_species_length_summary.csv"))
wigley_sizes_shelf <- read_csv(here::here("Data/processed/shelfwide_wigley_species_size_summary.csv"))




# Join them
size_results <- bind_rows(
  list(
    "Finfish Community" = bind_rows(finfish_size_df, finfish_lengths_shelf),
    "Wigley Subset" = bind_rows(wigley_size_df, wigley_sizes_shelf)), 
  .id = "community"
) %>% 
  mutate(
    survey_area = factor(survey_area, levels = area_levels_all)
  )




# Get trends:
len_mod_ffish <- lmer(
  formula = med_len_cm ~ scale(est_year) * survey_area * season + (1|est_year),
  data = finfish_size_df)
# len_mod_wig <- lmer(
#   formula = med_len_cm ~ est_year * survey_area * season + (1|est_year),
#   data = wigley_size_df)
wt_mod_wig <- lmer(
  formula = med_wt_kg ~ scale(est_year) * survey_area * season + (1|est_year),
  data = drop_na(wigley_size_df, med_wt_kg))



# Get the predictions and flag whether they are significant
size_trend_predictions <- bind_rows(list(
  "med_len_cm-Finfish Community" = get_preds_and_trendsignif(len_mod_ffish),
  "med_wt_kg-Wigley Subset"   = get_preds_and_trendsignif(wt_mod_wig)
  ), .id = "model_id") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric))






# Left side: 
# Median length - all species
# color = season vs. annual
# facet_wrap(~survey_area)
size_long <- size_results %>% 
  pivot_longer(cols = c(med_len_cm, med_wt_kg), 
               names_to = "metric", 
               values_to = "val")  %>% 
  mutate(survey_area = fct_relevel(survey_area, area_levels_all),
         season = fct_rev(season))


# just length for scales
len_plot <-  size_long %>% 
  filter(metric == "med_len_cm",
         community == "Finfish Community") %>% 
  drop_na(val) %>% 
  ggplot() +
  geom_ribbon(
    data = filter(size_trend_predictions, metric == "med_len_cm", non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    #data = filter(size_long, metric == "med_len_cm"),
    aes(est_year, val, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, val, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  geom_line(
    data = filter(size_trend_predictions, metric == "med_len_cm", non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(survey_area ~ metric*community,
             scales = "free") +
  theme(strip.text.y = element_blank()) +
  labs(
    x = "Year",
    y = "Length (cm)",
    color = "Season",
    linetype = "Linetype",
    fill = "Season")

# Right: 
# Median weight - wigley species
# color = season vs. annual
# facet_wrap(~survey_area)

# weight plots
wt_plot <- size_long %>% 
  filter(metric == "med_wt_kg") %>% 
  drop_na(val) %>% 
  ggplot() +
  geom_ribbon(
    data = filter(size_trend_predictions, metric == "med_wt_kg", non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    aes(est_year, val, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, val, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  geom_line(
    data = filter(size_trend_predictions, metric == "med_wt_kg", non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(survey_area ~ metric*community,
             scales = "free")  +
  labs(
    x = "Year",
    y = "Weight (kg)",
    color = "Season",
    linetype = "Linetype",
    fill = "Season")

(len_plot | wt_plot) + plot_layout(guides = "collect") 

# (len_plot | wt_plot) + plot_layout(guides = "collect", widths = c(3,1.5)) & theme(legend.position = "bottom")

Notes:

There are a number of years when median weight surges in MAB. For these situations I think we likely could trace these down to exceptionally high catches of larger species like elasmobranchs.

We also don’t see a shelf-wide decline in size so I might need to check what Friedland did or prepare for a confrontation.

C. Small/Large Fish Trends

The following figures frame “large fish” as individuals with length/bodymass greater than the 90th percentile for all years within that region. “Small fish” are considered smaller than the 10th percentile of the distribution.

Originally this had been done independently of species, I have now changed that

NOTE: Percentiles are determined using DescTools::Quantile() which uses numlen_adj to weight everything by abundance.

All Finfish Length LFI

If we look at the size distribution across all species these are the cutoff percentiles based on length.

# Get percentiles based on length for all species
DescTools::Quantile(
  finfish_in$length_cm, 
  weights = finfish_in$numlen_adj, 
  probs = c(0.05, 0.1, 0.5, 0.90, .95)) %>% 
      t() %>% 
      as_tibble() %>% 
  mutate(survey_area = "Northeast Shelf", .before = "5%") %>% 
  gt() %>% 
  tab_header(
    title = "All Finfish Species - Body Length Quantiles (cm)") %>% 
  tab_spanner(label = "Percentile", columns = -survey_area) %>% 
  cols_label("survey_area" = "")

All Finfish Species - Body Length Quantiles (cm)
	Percentile
	5%	10%	50%	90%	95%
Northeast Shelf	4	5	14	41	66

Looking instead at the size distributions independently for all species these are the new cutoff percentiles based on length for ten species.

# Get percentiles based on length for each species
# # Get 90th or 95th percentile as size threshold
# Do them separately
species_size_quants <- finfish_in %>% 
  mutate(survey_area = "Northeast Shelf") %>% 
  unite("species_area", c(survey_area, comname), sep = "X", ) %>% 
  split(.$species_area) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      x$length_cm, 
      #x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.05, 0.1, 0.5, 0.90, .95)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "species_area") %>% 
  separate(col = species_area, into = c("survey_area", "comname"), sep = "X") %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels_all)) %>% 
  arrange(survey_area)



# Show the table
species_size_quants %>% 
  mutate_if(is.numeric, round, 2) %>% 
  head(n = 10) %>% 
  gt() %>% 
  tab_header(title = "Finfish Community - Body Length Quantiles (cm)") %>% tab_spanner(label = "Percentile", columns = -survey_area) %>% 
  cols_label("survey_area" = "")

Finfish Community - Body Length Quantiles (cm)
	Percentile
	comname	5%	10%	50%	90%	95%
Northeast Shelf	acadian redfish	14.00	17.0	25.0	32.0	34.0
Northeast Shelf	african pompano	4.20	4.4	6.0	7.6	7.8
Northeast Shelf	albacore	66.00	66.0	66.0	66.0	66.0
Northeast Shelf	alewife	11.00	12.0	20.0	25.0	27.0
Northeast Shelf	alligatorfish	7.00	7.0	11.0	13.0	14.0
Northeast Shelf	almaco jack	26.00	26.0	26.0	26.0	26.0
Northeast Shelf	american eel	29.05	30.1	35.5	53.8	59.7
Northeast Shelf	american plaice	13.00	15.0	26.0	39.0	44.0
Northeast Shelf	american sand lance	14.00	14.0	16.0	18.0	19.0
Northeast Shelf	american shad	17.00	19.0	24.0	33.0	38.0

The following figure shows changes in seasonal large and small fish indices within each region. The large and small size thresholds are set for each species independently using their length distributions as in the table above.

The axes on these are: \[\frac{\sum{individuals_{ijsk} > size~threshold_k}} {\sum{individuals_{ijs}}} * 100\]

for each: \(year_i\), \(region_j\), \(season_s\) & \(species_k\).

# Assign large/small indicators
finfish_lfi <- finfish_in %>% 
  left_join(select(species_size_quants, comname, `10%`, `90%`)) %>%   
  mutate(
    small = ifelse(length_cm < `10%`, T, F),
    large = ifelse(length_cm > `90%`, T, F))


# Get total abundance in each area*year*season
total_abundbio <- bind_rows(
    mutate(finfish_lfi, survey_area = "Northeast Shelf"),
    finfish_lfi) %>% 
  group_by(survey_area, est_year, season) %>% 
  summarise(
    total_abund = sum(numlen_adj, na.rm = T),
    .groups = "drop")
  
  
# Get large/small fish  abundances
lf_abundbio <- bind_rows(
    mutate(finfish_lfi, survey_area = "Northeast Shelf"),
    finfish_lfi) %>%
  group_by(survey_area, comname, est_year, season) %>%
  summarise(
    small_abund = sum(numlen_adj * small),
    large_abund = sum(numlen_adj * large),
    .groups = "drop") %>% 
  group_by(survey_area, est_year, season) %>% 
  summarise(
    small_abund = sum(small_abund, na.rm = T),
    large_abund = sum(large_abund, na.rm = T),
    .groups = "drop")

# Join and divide
lfi <- left_join(total_abundbio, lf_abundbio) %>%
  mutate(
    LFI = large_abund / total_abund,
    SFI = small_abund / total_abund)  %>%
  mutate(
    survey_area = factor(survey_area, levels = area_levels_all),
    season = factor(season, levels = c("Spring", "Fall"))) 


# #plot it
# lfi %>%
#   pivot_longer(cols = c(LFI, SFI), names_to = "Index", values_to = "val") %>%
#   mutate(size_thresh = if_else(Index == "LFI", ">90th Percentile\nSpecies-Based", "<10th Percentile\nSpecies-Based")) %>%
#   ggplot(aes(est_year, val, color = season)) +
#   geom_point(
#     aes(est_year, val, color = season),
#     size = 0.8, alpha = 0.5) +
#   geom_ma(
#     aes(est_year, val, color = season, linetype = "5-Year Moving Average"), 
#     n = 5, ma_fun = SMA) +
#   scale_fill_gmri() +
#   scale_color_gmri() +
#   scale_y_continuous(labels = label_percent()) +
#   facet_nested(survey_area~fct_rev(Index)*size_thresh) +
#   labs(
#     x = "Year",
#     y = "Percent of Overall Abundance",
#     title = "All Finfish, Large and Small Fish Indices, by Length")

# Plot them as loing-term average

# Plot the SFI
sfi_p <- lfi %>% 
  ggplot() +
  geom_point(
    aes(est_year, SFI, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, SFI, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(labels = label_percent()) +
  theme(strip.text.y = element_blank()) +
  facet_grid(survey_area~., scales = "free") +
  labs(
    x = "Year",
    y = "Percent of Total Abundance",
    title = "Finfish Species",
    subtitle = "Small Fish Index\n(<10 Percentile\nSpecies-Length Distribution)",
    fill = "Season",
    color = "Season",
    linetype = "")


# Plot the LFI
lfi_p <- ggplot(lfi) +
  geom_point(
    aes(est_year, LFI, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, LFI, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(labels = label_percent()) +
  facet_grid(survey_area~., scales = "free") +
  labs(
    x = "Year",
    y = "",
    subtitle = "Large Fish Index\n(>90 Percentile\nSpecies-Length Distribution)",
    fill = "Season",
    color = "Season",
    linetype = "")



# Plot together
(sfi_p | lfi_p) + plot_layout(guides = "collect")

Wigley Regional LFI

There was a question on whether to do LFI and SFI based on the species quantiles, or simply overall regional quantiles.

The following table shows what the bodymass percentiles are for each region for the Wigley subset. These ignore seasonal and species specific benchmarks.

# # Get 90th or 95th percentile as size threshold
# Do them separately
area_size_quants <- bind_rows(
    mutate(wigley_in, survey_area = "Northeast Shelf"),
    wigley_in) %>% 
  split(.$survey_area) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      # x$length_cm, 
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.05, 0.1, 0.5, 0.90, .95)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "survey_area") %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels_all)) %>% 
  arrange(survey_area)

# Show the table
area_size_quants %>% 
  mutate_if(is.numeric, round, 2) %>% 
  gt() %>% 
  tab_header(title = "Wigley Species - Body Mass Quantiles (kg)") %>% tab_spanner(label = "Percentile", columns = -survey_area) %>% 
  cols_label("survey_area" = "")

Wigley Species - Body Mass Quantiles (kg)
	Percentile
	5%	10%	50%	90%	95%
Northeast Shelf	0.01	0.01	0.08	0.93	1.61
GoM	0.01	0.02	0.14	1.11	1.68
GB	0.00	0.01	0.14	0.91	1.51
SNE	0.00	0.01	0.05	0.58	1.43
MAB	0.01	0.02	0.04	1.21	1.68

Then this figure shows what the LFI and SFI indices are using those same overall benchmarks which ignore species:

# Assign large/small
wigley_in <- wigley_in %>% 
  left_join(select(area_size_quants, survey_area, `10%`, `90%`)) %>%   
  mutate(
    small = ifelse(ind_weight_kg < `10%`, T, F),
    large = ifelse(ind_weight_kg > `90%`, T, F))



# Get total abundance
total_abundbio <- bind_rows(
    mutate(wigley_in, survey_area = "Northeast Shelf"),
    wigley_in) %>% 
  group_by(survey_area, est_year, season) %>% 
  summarise(
    total_abund = sum(numlen_adj),
    total_bio = sum(sum_weight_kg),
    .groups = "drop")
  
  
# Get large fish bio
lf_abundbio <- bind_rows(
    mutate(wigley_in, survey_area = "Northeast Shelf"),
    wigley_in) %>% 
  group_by(survey_area, est_year, season) %>% 
  summarise(
    small_abund = sum(numlen_adj * small),
    large_abund = sum(numlen_adj * large),
    .groups = "drop")



# Join and divide
lfi <- left_join(total_abundbio, lf_abundbio) %>% 
  mutate(
    LFI = large_abund / total_abund,
    SFI = small_abund / total_abund)  %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels_all)) %>% 
  left_join(
    select(area_size_quants, survey_area, `10%`, `90%`) %>% 
    mutate(across(-survey_area, \(x) str_c(round(x, 2), " kg")))) %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall")))

# Plot them as loing-term average
# Plot the SFI
sfi_p <- lfi %>% 
  ggplot() +
  geom_point(
    aes(est_year, SFI, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, SFI, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(labels = label_percent()) +
  theme(strip.text.y = element_blank()) +
  facet_grid(survey_area~., scales = "free") +
  labs(
    x = "Year",
    y = "Percent of Total Abundance",
    title = "Wigley Species",
    subtitle = "Small Fish Index\n(<10 Percentile\nRegional-Mass Distribution)",
    fill = "Season",
    color = "Season",
    linetype = "")


# Plot the LFI
lfi_p <- ggplot(lfi) +
  geom_point(
    aes(est_year, LFI, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, LFI, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(labels = label_percent()) +
  facet_grid(survey_area~., scales = "free") +
  labs(
    x = "Year",
    y = "",
    subtitle = "Large Fish Index\n(>90 Percentile\nRegional-Mass Distribution)",
    fill = "Season",
    color = "Season",
    linetype = "")



# Plot together
(sfi_p | lfi_p) + plot_layout(guides = "collect")

Using these broad group quantiles, we have the following species are represented within the large fish index:

big_fishes <- wigley_in %>% 
  filter(large) %>% 
  distinct(comname) %>% 
  arrange(comname) %>% 
  rename(`Common Name:` = comname) 


# Print it neatly
pander(big_fishes$`Common Name:`)

acadian redfish, alewife, american plaice, american shad, atlantic angel shark, atlantic cod, atlantic croaker, atlantic halibut, atlantic mackerel, atlantic sharpnose shark, atlantic sturgeon, atlantic torpedo, atlantic wolffish, barndoor skate, black sea bass, blackbelly rosefish, blueback herring, bluefish, bluntnose stingray, buckler dory, bullnose ray, butterfish, clearnose skate, cownose ray, cunner, cusk, fawn cusk-eel, fourspot flounder, goosefish, greater amberjack, haddock, little skate, longhorn sculpin, northern searobin, ocean pout, offshore hake, pollock, red hake, roughtail stingray, sand tiger, sandbar shark, scup, sea raven, silver hake, smooth butterfly ray, smooth dogfish, smooth skate, southern stingray, spanish mackerel, spiny butterfly ray, spiny dogfish, spotted hake, striped bass, summer flounder, tautog, thorny skate, weakfish, white hake, windowpane flounder, winter flounder, winter skate, witch flounder and yellowtail flounder

You might think: “Why is a 1kg acadian redfish a”large fish”?“. That’s a consequence of abundance declining with size pretty dramatically.

There’s really just so few “large” individuals that the percentiles sit quite small.

Wigley Species-Specific LFI

Using species-specific size distribution quantiles shifts the story very little: so little that I feel like its bugged, but I can’t track down where

# Get percentiles based on length for each species
# # Get 90th or 95th percentile as size threshold
# Do them separately
species_size_quants <- wigley_in %>% 
  mutate(survey_area = "Northeast Shelf") %>% 
  unite("species_area", c(survey_area, comname), sep = "X", ) %>% 
  split(.$species_area) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      # x$length_cm, 
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.05, 0.1, 0.5, 0.90, .95)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "species_area") %>% 
  separate(col = species_area, into = c("survey_area", "comname"), sep = "X") %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels_all)) %>% 
  arrange(survey_area)


# Show the table
species_size_quants %>% 
  mutate_if(is.numeric, round, 2) %>% 
  head(n = 10) %>% 
  gt() %>% 
  tab_header(title = "Wigley Community - Body Mass Quantiles (kg)") %>% tab_spanner(label = "Percentile", columns = -survey_area) %>% 
  cols_label("survey_area" = "")

Wigley Community - Body Mass Quantiles (kg)
	Percentile
	comname	5%	10%	50%	90%	95%
Northeast Shelf	acadian redfish	0.04	0.07	0.25	0.55	0.67
Northeast Shelf	alewife	0.01	0.01	0.09	0.21	0.28
Northeast Shelf	american plaice	0.01	0.02	0.14	0.53	0.76
Northeast Shelf	american shad	0.05	0.07	0.15	0.42	0.66
Northeast Shelf	atlantic angel shark	0.33	0.43	1.45	9.34	14.92
Northeast Shelf	atlantic cod	0.04	0.15	1.21	4.92	6.84
Northeast Shelf	atlantic croaker	0.08	0.10	0.15	0.34	0.41
Northeast Shelf	atlantic halibut	0.28	0.40	1.66	8.60	13.71
Northeast Shelf	atlantic herring	0.01	0.02	0.09	0.16	0.17
Northeast Shelf	atlantic mackerel	0.05	0.05	0.12	0.30	0.37

The following figure shows changes in seasonal large and small fish indices within each region. The large and small size thresholds are set for each species independently using their weight distributions as in the table above.

# Assign large/small indicators
wigley_lfi <- wigley_in %>% 
  select(-small, -large) %>% 
  left_join(select(species_size_quants, comname, `10%`, `90%`)) %>%   
  mutate(
    small = ifelse(ind_weight_kg < `10%`, T, F),
    large = ifelse(ind_weight_kg > `90%`, T, F))


# Get total abundance in each area*year*season
total_abundbio <- bind_rows(
    mutate(wigley_lfi, survey_area = "Northeast Shelf"),
    wigley_lfi) %>% 
  group_by(survey_area, est_year, season) %>% 
  summarise(
    total_abund = sum(numlen_adj, na.rm = T),
    .groups = "drop")
  
  
# Get large/small fish  abundances
lf_abundbio <- bind_rows(
    mutate(wigley_lfi, survey_area = "Northeast Shelf"),
    wigley_lfi) %>%
  group_by(survey_area, comname, est_year, season) %>%
  summarise(
    small_abund = sum(numlen_adj * small),
    large_abund = sum(numlen_adj * large),
    .groups = "drop") %>% 
  group_by(survey_area, est_year, season) %>% 
  summarise(
    small_abund = sum(small_abund, na.rm = T),
    large_abund = sum(large_abund, na.rm = T),
    .groups = "drop")

# Join and divide
lfi <- left_join(total_abundbio, lf_abundbio) %>%
  mutate(
    LFI = large_abund / total_abund,
    SFI = small_abund / total_abund)  %>%
  mutate(
    survey_area = factor(survey_area, levels = area_levels_all),
    season = factor(season, levels = c("Spring", "Fall"))) 

# Plot them as loing-term average

# Plot the SFI
sfi_p <- lfi %>% 
  ggplot() +
  geom_point(
    aes(est_year, SFI, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, SFI, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(labels = label_percent()) +
  theme(strip.text.y = element_blank()) +
  facet_grid(survey_area~., scales = "free") +
  labs(
    x = "Year",
    y = "Percent of Total Abundance",
    title = "Wigley Species",
    subtitle = "Small Fish Index\n(<10 Percentile\nSpecies-Mass Distribution)",
    fill = "Season",
    color = "Season",
    linetype = "")


# Plot the LFI
lfi_p <- ggplot(lfi) +
  geom_point(
    aes(est_year, LFI, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, LFI, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(labels = label_percent()) +
  facet_grid(survey_area~., scales = "free") +
  labs(
    x = "Year",
    y = "",
    subtitle = "Large Fish Index\n(>90 Percentile\nSpecies-Mass Distribution)",
    fill = "Season",
    color = "Season",
    linetype = "")



# Plot together
(sfi_p | lfi_p) + plot_layout(guides = "collect")

Groundfish Complex LFI

For this section I’m going to do something a little different to dial in to species specific body size changes.

For the groundfish complex species, across all individuals estimate weight percentiles, and use those to get a LFI and SFI for those species

# # Get 90th or 95th percentile as size threshold
# Do them separately
species_size_quants <- gfish_comp %>% 
  mutate(survey_area = "Northeast Shelf") %>% 
  unite("species_area", c(survey_area, comname), sep = "X", ) %>% 
  split(.$species_area) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      # x$length_cm, 
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.05, 0.1, 0.5, 0.90, .95)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "species_area") %>% 
  separate(col = species_area, into = c("survey_area", "comname"), sep = "X") %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels_all)) %>% 
  arrange(survey_area)



# # Show the table
# species_size_quants %>% 
#   mutate_if(is.numeric, round, 2) %>% 
#   gt() %>% 
#   tab_header(title = "Groundfish Complex - Body Mass Quantiles (kg)") %>% tab_spanner(label = "Percentile", columns = -survey_area) %>% 
#   cols_label("survey_area" = "")



# Assign large/small
gfish_comp <- gfish_comp %>% 
  left_join(select(species_size_quants, comname, `10%`, `90%`)) %>%   
  mutate(
    small = ifelse(ind_weight_kg < `10%`, T, F),
    large = ifelse(ind_weight_kg > `90%`, T, F))



# Get total abundance
total_abundbio <- bind_rows(
    mutate(gfish_comp, survey_area = "Northeast Shelf"),
    gfish_comp) %>% 
  group_by(survey_area, est_year, season) %>% 
  summarise(
    total_abund = sum(numlen_adj),
    total_bio = sum(sum_weight_kg),
    .groups = "drop")
  
  
# Get large fish bio
lf_abundbio <- bind_rows(
    mutate(gfish_comp, survey_area = "Northeast Shelf"),
    gfish_comp) %>% 
  group_by(survey_area, est_year, season) %>% 
  summarise(
    small_abund = sum(numlen_adj * small),
    small_bio = sum(sum_weight_kg * small),
    large_abund = sum(numlen_adj * large),
    large_bio = sum(sum_weight_kg * large),
    .groups = "drop")

# Join and divide
lfi <- left_join(total_abundbio, lf_abundbio) %>% 
  mutate(
    LFI = large_abund / total_abund,
    SFI = small_abund / total_abund)  %>% 
  mutate(
    survey_area = factor(survey_area, levels = area_levels_all),
    season = factor(season, levels = c("Spring", "Fall"))) %>% 
  left_join(
    select(area_size_quants, survey_area, `10%`, `90%`) %>% 
      mutate(across(-survey_area, \(x) str_c(round(x, 2), " kg"))))

# Plot them as loing-term average

# Plot the SFI
sfi_p <- lfi %>% 
  ggplot() +
  geom_point(
    aes(est_year, SFI, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, SFI, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(labels = label_percent()) +
  theme(strip.text.y = element_blank()) +
  facet_grid(survey_area~., scales = "free") +
  labs(
    x = "Year",
    y = "Percent of Total Abundance",
    title = "Groundfish Complex",
    subtitle = "Small Fish Index\n(<10 Percentile\nSpecies-Mass Distribution)",
    fill = "Season",
    color = "Season",
    linetype = "")


# Plot the LFI
lfi_p <- ggplot(lfi) +
  geom_point(
    aes(est_year, LFI, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, LFI, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(labels = label_percent()) +
  facet_grid(survey_area~., scales = "free") +
  labs(
    x = "Year",
    y = "",
    subtitle = "Large Fish Index\n(>90 Percentile\nSpecies-Mass Distribution)",
    fill = "Season",
    color = "Season",
    linetype = "")



# Plot together
(sfi_p | lfi_p) + plot_layout(guides = "collect")

Shift in Regional Percentile Size

A 90th percentile fish is smaller than I feel it has been in the past, so I wanted to plot those quantiles over time.

Indeed, the 90th & 95th percentile size for fishes in the GOM was nearly 2x what it is today.

# Do them separately
yearly_size_quants <- bind_rows(
    mutate(wigley_in, survey_area = "Northeast Shelf"),
    wigley_in) %>% 
  unite(col = "group_var", survey_area, est_year, sep = "X", remove = T) %>% 
  split(.$group_var) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      # x$length_cm, 
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.90, .95)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "group_var") %>% 
  separate(group_var, into = c("survey_area", "est_year"), sep = "X") %>% 
  mutate(
    est_year = as.numeric(est_year),
    survey_area = factor(survey_area, levels = area_levels_all)) %>% 
  arrange(survey_area)


# Reshape to make filtering easier
quants_long <- yearly_size_quants %>% 
  #select(-c(`5%`,`50%`, `95%`)) %>% 
  pivot_longer(cols = -c(survey_area, est_year), names_to = "Percentile", values_to = "biomass_kg") 


# Separate axes
q10_p <- quants_long %>%
  filter(Percentile == "10%") %>% 
  ggplot(aes(est_year, biomass_kg)) +
  geom_point(size = 0.8, alpha = 0.5) +
  geom_ma(n = 5, ma_fun = SMA, aes(linetype = "5-Year Moving Average"), color = "black") +
  scale_color_gmri() +
  theme(
    strip.text.y = element_blank(),
    plot.margin = margin(0,0,0,0)) +
  facet_grid(survey_area~Percentile) +
  labs(y = "Bodymass (kg)", x = "Year", linetype = "")

# Separate axes
q90_p <- quants_long %>%
  filter(Percentile == "90%") %>% 
  ggplot(aes(est_year, biomass_kg)) +
  geom_point(size = 0.8, alpha = 0.5) +
  geom_ma(n = 5, ma_fun = SMA, aes(linetype = "5-Year Moving Average"), color = "black") +
  scale_color_gmri(reverse = T) +
  theme(
    strip.text.y = element_blank(),
    plot.margin = margin(0,0,0,0)) +
  facet_grid(survey_area~Percentile) +
  labs(y = "", x = "Year", linetype = "")


# Separate axes
q95_p <- quants_long %>%
  filter(Percentile == "95%") %>% 
  ggplot(aes(est_year, biomass_kg)) +
  geom_point(size = 0.8, alpha = 0.5) +
  geom_ma(n = 5, ma_fun = SMA, aes(linetype = "5-Year Moving Average"), color = "black") +
  scale_color_gmri(reverse = T) +
  theme(plot.margin = margin(0,0,0,0)) +
  facet_grid(survey_area~Percentile) +
  labs(y = "", x = "Year", linetype = "")

(q10_p | q90_p | q95_p) + plot_layout(guides = "collect", axes = "collect_x") +
  plot_annotation(title = "Size Percentile Changes")

10th percentile weights seem stable. 90th percentile body-weights have fluctuated over time.

D. Decadal Binned Spectra Plot

The following figure is a visual aid of the normalized biomass spectra for each decade.

Spring

# We need to normalize to be comparable to the Edwards method
# Might be able to just assign the bins and their widths, then divide the weights by that
# ^ seems to work!

# add the spectra slope values from MLE method as a side value, could use facet grid label

# plot the normalized
wigley_in %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels)) %>% 
  filter(ind_weight_kg <= 30,
         season == "Spring") %>% 
  assign_log2_bins(log2_increment = 1) %>% 
  mutate(season = season,
         decade = str_c(floor_decade(est_year), "'s")) %>% 
  ggplot(aes(
    x = ind_weight_g, weight = numlen_adj/bin_width)) + 
  geom_histogram(
    aes(fill = after_stat(density)),
    breaks = 2^seq(0,16,1), color = "gray90", linewidth = 0.2) +
  scale_y_continuous(
    transform = transform_log(base = 2),
    labels = label_log(base = 2),
    expand = expansion(add = c(0,0))) +
  scale_x_continuous(
    transform = transform_log(base = 2),
    labels = label_log(base = 2),
    breaks = 2^seq(0,16, 2),
    limits = c(1, 2^14),
    expand = expansion(add = c(0,0))) +
  scale_fill_distiller(palette = "YlGnBu", direction = 1) +
  scale_color_gmri() +
  #facet_nested(decade~survey_area*fct_rev(season)) +
  facet_nested(survey_area*fct_rev(season)~decade) +
  theme(
    axis.text.y = element_blank(), 
    panel.grid.major.y = element_blank(), 
    legend.position = "bottom") +
  labs(
    title = "Normalized Biomass Spectra",
    x = "Mass (g)",
    y = "Density",
    fill = "Normalized Abundance Density\n(abundance / bin-width)")

Fall

# We need to normalize to be comparable to the Edwards method
# Might be able to just assign the bins and their widths, then divide the weights by that
# ^ seems to work!

# add the spectra slope values from MLE method as a side value, could use facet grid label

# plot the normalized
wigley_in %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels)) %>% 
  filter(ind_weight_kg <= 30,
         season == "Fall") %>% 
  assign_log2_bins(log2_increment = 1) %>% 
  mutate(season = season,
         decade = str_c(floor_decade(est_year), "'s")) %>% 
  ggplot(aes(
    x = ind_weight_g, weight = numlen_adj/bin_width)) + 
  geom_histogram(
    aes(fill = after_stat(density)),
    breaks = 2^seq(0,16,1), color = "gray90", linewidth = 0.2) +
  scale_y_continuous(
    transform = transform_log(base = 2),
    labels = label_log(base = 2),
    expand = expansion(add = c(0,0))) +
  scale_x_continuous(
    transform = transform_log(base = 2),
    labels = label_log(base = 2),
    breaks = 2^seq(0,16, 2),
    limits = c(1, 2^14),
    expand = expansion(add = c(0,0))) +
  scale_fill_distiller(palette = "YlGnBu", direction = 1) +
  scale_color_gmri() +
  #facet_nested(decade~survey_area*fct_rev(season)) +
  facet_nested(survey_area*fct_rev(season)~decade) +
  theme(
    axis.text.y = element_blank(), 
    panel.grid.major.y = element_blank(), 
    legend.position = "bottom") +
  labs(
    title = "Normalized Biomass Spectra",
    x = "Mass (g)",
    y = "Density",
    fill = "Normalized Abundance Density\n(abundance / bin-width)")

Spring & Fall

# We need to normalize to be comparable to the Edwards method
# Might be able to just assign the bins and their widths, then divide the weights by that
# ^ seems to work!

# add the spectra slope values from MLE method as a side value, could use facet grid label

# plot the normalized
wigley_in %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels)) %>% 
  filter(ind_weight_kg <= 30) %>% 
  assign_log2_bins(log2_increment = 1) %>% 
  mutate(season = season,
         decade = str_c(floor_decade(est_year), "'s")) %>% 
  ggplot(aes(
    x = ind_weight_g, weight = numlen_adj/bin_width)) + 
  geom_histogram(
    aes(fill = after_stat(density)),
    breaks = 2^seq(0,16,1), color = "gray90", linewidth = 0.2) +
  scale_y_continuous(
    transform = transform_log(base = 2),
    labels = label_log(base = 2),
    expand = expansion(add = c(0,0))) +
  scale_x_continuous(
    transform = transform_log(base = 2),
    labels = label_log(base = 2),
    breaks = 2^seq(0,16, 2),
    limits = c(1, 2^14),
    expand = expansion(add = c(0,0))) +
  scale_fill_distiller(palette = "YlGnBu", direction = 1) +
  scale_color_gmri() +
  #facet_nested(decade~survey_area*fct_rev(season)) +
  facet_nested(survey_area*fct_rev(season)~decade) +
  theme(
    axis.text.y = element_blank(), 
    panel.grid.major.y = element_blank(), 
    legend.position = "bottom") +
  labs(
    title = "Normalized Biomass Spectra",
    x = "Mass (g)",
    y = "Density",
    fill = "Normalized Abundance Density\n(abundance / bin-width)")

Summary of Trends

Let’s start with these pieces, and we can then discuss doing something with the extremes of the size distribution–changes in small vs. large individuals, where we pre-define size bins as you’ve done in slides 7 and 8, or perhaps using an approach similar to Lora’s approach of defining the upper/lower percentiles of size. - KMillz

---

Part 2: Relationships to External Factors

The aim of this section is to provide context and potentially attribution/correlations that may help explain trends/changes observed in the above section.

I will lead with a characterization of broad regional changes in temperature and fisheries landings totals.

# # Model Dataset for Wigley Species Subset
wigley_bmspectra_model_df <- read_csv(here::here("Data/model_ready/wigley_community_bmspectra_mod.csv"))

# load shelfwide stuff too
wigley_bmspectra_shelf <- read_csv(here::here("Data/processed/shelfwide_wigley_species_bodymass_spectra.csv"))
bot_temps <- read_csv(here::here("Data/processed", "trawl_region_seasonal_bottom_temps.csv"))


# Bolt that on
wigley_bmspectra_model_df <- bind_rows(
  wigley_bmspectra_model_df,
  left_join(wigley_bmspectra_shelf, bot_temps, by  = join_by("est_year" == "year", survey_area, season)))


# Use Wigley bmass spectra model data for plots
wigley_bmspectra_model_df <- wigley_bmspectra_model_df %>% 
  mutate(survey_area = case_when(
    survey_area == "GoM" ~ "Gulf of Maine",
    survey_area == "GB" ~ "Georges Bank",
    survey_area == "SNE" ~ "Southern New England",
    survey_area == "MAB" ~ "Mid-Atlantic Bight",
    survey_area == "Northeast Shelf" ~"Northeast Shelf"),
  survey_area = factor(survey_area, levels = area_levels_long_all),
  season = fct_rev(season))

Bottom Temperature

Bottom Temperature Trends

Bottom temperature data is from the DuPontavice and Saba combined bottom temperature product. This data product has been clipped to the study areas and seasonally averaged to match the survey sampling seasons.

Plotted below are the long-term trends in each area for each survey season. Significant seasonal trends are shown with linear trends and error displayed.

# Model linear trends in time
btemp_mods <- lmer(
  formula = bot_temp ~ scale(est_year) * survey_area * season + (1|est_year),
  data = wigley_bmspectra_model_df)



# change the stupid labels around because
btemp_preds <- as.data.frame(
    # Model predictions
    ggpredict(
        btemp_mods, 
        terms = ~ est_year * survey_area * season) ) %>% 
    mutate(
        survey_area = factor(group, levels = area_levels_long_all),
        season = factor(facet, levels = c("Spring", "Fall")))
btemp_signif <- emtrends(
      object = btemp_mods,
      specs =  ~ survey_area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
btemp_trend_predictions <- left_join(btemp_preds, btemp_signif)  %>% 
  mutate(survey_area = group,
         survey_area = factor(survey_area, levels = area_levels_long_all),
         season = factor(season, levels = c("Spring", "Fall")))




# Plot temperatures
temp1 <- wigley_bmspectra_model_df %>% 
  ggplot() +
  geom_ribbon(
    data = filter(btemp_trend_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    aes(est_year, bot_temp, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, bot_temp, color = season, linetype = "5-Year Smooth"),
    size = 1, n = 5, ma_fun = SMA) +
  geom_line(
    data = filter(btemp_trend_predictions, non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_color_gmri() + 
  scale_fill_gmri() + 
  #scale_linetype_manual(values = c(1,1)) +
  facet_grid(survey_area~., scales = "free") +
  theme(legend.position = "right") +
  scale_y_continuous(labels = label_number(suffix = "\u00b0C")) +
  labs(y = "Bottom Temperature", x = "Year", color = "Season",
       fill = "Season", linetype = "Trend",
       title = "Significant trends in Regional Bottom Temperature")



temp1 

Seasonal bottom temperatures have increased in both spring and fall for the Northeast Shelf.

Looking at different areas within the region, Fall temperatures are increasing across the board, but seasonally spring temperatures in SNE and MAB are stable.

Seasonal Bottom Temperature Ranges

The seasonal difference in bottom temperatures are such that the more southern regions experience a larger spring-fall difference in absolute temperatures each year.

Within the Gulf of Maine there is a notably small difference when compared against the other regions. The seasonal difference in bottom temperatures here can be smaller than the North-to-South gradient in spring bottom temperatures across the whole Northeast Shelf.

The overlap in bottom temperature distributions is very similar to the overlap in mass spectra exponents. Gulf of Maine spectra and bottom temperatures share a good region of overlap, and the regions that experience larger seasonal fluctuations in BT see larger seasonal fluctuations in their size distributions.

# Plot the distribution as a raincloud plot
seasonal_temp_dist <- wigley_bmspectra_model_df %>% 
  ggplot(aes(fct_rev(survey_area), bot_temp, color = season, fill = season)) +
   ggdist::stat_halfeye(
     alpha = 0.4,
     adjust = .5, 
     width = .6, 
    .width = 0, 
    justification = -.2, 
    point_colour = NA) + 
  geom_boxplot(
    fill = "transparent",
    width = .15, 
    outlier.shape = NA) +
  ## add dot plots from {ggdist} package
  ggdist::stat_dots(
    ## orientation to the left
    side = "left", size = 1,
    ## move geom to the left
    justification = 1.12, 
    ## adjust grouping (binning) of observations 
    binwidth = .1) + 
  scale_color_gmri() + 
  scale_fill_gmri() + 
  scale_y_continuous(labels = label_number(suffix = "\u00b0C")) +
  coord_flip() +
  theme(plot.margin = margin(c(0,0,0,0))) +
  theme(legend.position = "none") +
  labs(y = "Bottom Temperature", x = "", color = "Season",
       title = "Seasonal Bot Temperature")


# Distribution of size distribution exponenents
wigley_b_dist_plot <- wigley_bmspectra_model_df %>% 
  ggplot(aes(fct_rev(survey_area), b, color = season, fill = season)) +
  ggdist::stat_halfeye(
     alpha = 0.4,
     adjust = .5, 
     width = .6, 
    .width = 0, 
    justification = -.2, 
    point_colour = NA) + 
  geom_boxplot(
    fill = "transparent",
    width = .15, 
    outlier.shape = NA) +
  ## add dot plots from {ggdist} package
  ggdist::stat_dots(
    ## orientation to the left
    side = "left", size = 1,
    ## move geom to the left
    justification = 1.12, 
    ## adjust grouping (binning) of observations 
    binwidth = .005) + 
  scale_color_gmri() + 
  scale_fill_gmri() + 
  coord_flip() +
  theme(plot.margin = margin(c(0,0,0,0))) +
  labs(
    y = "Exponent of Size Distribution", 
    x = "",
    title = "Mass Distributions",
    subtitle = "Bodymass Ranges (min: 1g, max: max observed)")



# Plot next to bottom temperatures
(seasonal_temp_dist |  wigley_b_dist_plot) + plot_layout(guides = "collect")

Regional Temperature Relationship or Seasonal Catchability

The extent that the shape of the size distribution within each region follows some negative linear relationship to increasing temperature largely is a reflection of the seasonal differences.

If we follow the relationship within a specific season, the relationship is much weaker and not-significant in most cases.

wigley_bmspectra_model_df %>% 
  #filter(survey_area == "Northeast Shelf") %>% 
  ggplot(aes(bot_temp, b)) +
  geom_point(aes(color = season)) +
  geom_smooth(method = "lm", color = "black") +
  facet_wrap(~survey_area, ncol = 1) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(
    title = "Regional Temperature-Spectra Relationship\nReflects Seasonal Catchability Differences",
    subtitle = "Seasonal differences in spectra proportional to seasonal temperature range",
    shape = "Survey Area")

Temperature Sentiment

I believe that the size distribution along the shelf is responsive on shorter time scales to temperature changes than hypotheses that predict temperature related changes in growth. This is the seasonal variability in spectra we see.

This suggests to me that the community size distribution is likely less sensistive to long-term trends in temperature since the distribution is able to effectively track swings in temperature that happen within the year.

I think there is an important question about temperature’s role, but I don’t think that local warming is having much impact on the community size distribution.

If it is, it is secondary to other size distribution shifts happening at shorter time scales and across spatial scales.

This raises questions about how seasonal size distribution changes are being achieved.

From the decadal variability work we have reason to believe that individual species are both not tracking temperature changes in space at the rate they are ocurring AND that changes in size at age are happening. These are signals at the species level indicating that temperature effects on growth are ocurring.

There would need to be some compensatory mechanisms present to offset this across the broader community.

Framing “Warming” Locally via Local BT Anomalies

We need to avoid the conflation that a significant relationship between spectra and bottom temperature means that “warming” has altered the size spectra.

There are known relationships between inter-specific and intra-specific body size across temperature ranges and across latitudinal ranges.

• Bergmann’s Rule: Warmer climates are often inhabited by smaller-bodied species
• James Rule: At an intraspecific level, populations living in warmer conditions often reach smaller maximum body sizes

There is also metabolic theory that supports the basis for the above rules as well as seasonal redistribution related to body size and temperature.

• Metabolic Theory of Ecology: Metabolic theory predicts how metabolic rate, by setting the rates of resource uptake from the environment and resource allocation to survival, growth, and reproduction, controls ecological processes at all levels of organization from individuals to the biosphere.

Metabolic Efficiency and Seasonal Movement: MTE suggests that metabolic rates, which are influenced by body size and temperature, are central to the ecological behavior of organisms.

Larger individuals, which naturally have lower per-unit mass metabolic rates, may move to cooler areas seasonally to optimize their metabolic function and reduce the energetic costs of thermoregulation in warmer periods.

wigley_bmspectra_model_df %>% 
  filter(survey_area != "Northeast Shelf") %>% 
  ggplot(aes(bot_temp, b, color = season)) +
  geom_point(aes(shape = survey_area)) +
  geom_smooth(method = "lm", color = "black") +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  guides(shape = guide_legend(nrow = 2)) +
  theme(legend.box = "vertical") +
  labs(
    title = "Bergmann's Rule + MTE",
    subtitle = "Broad relationship between temperature and size distribution,\nstrengthened by seasonal temperature-related movements",
    shape = "Survey Area",
    color = "Season")

The observation that body size distribution broadly follows some temperature relationship is a confirmation of these broadly observed rules, but I would argue it does not convey that warming has shifted a local climate in some direction.

Since we know that individuals migrate and may use the full water column or occupy different habitats seasonally, we need to treat seasonal communities separately to avoid an apples to oranges comparison of the community size distribution.

Using these same season + region blocks, we can characterize the typical bottom temperature for these location and seasonal combinations using the long-term average, and quantify warming using departures from this long-term average.

This also resolves a separate issue when modeling with absolute temperatures, where temperature and season*area units are colinear.

Local Bottom Temperature Anomalies Model

This is the model I would suggest gets the most at “Warming” impacting regional size distributions:

Fixed effects:

• Season
• Survey area
• Bottom temperature anomaly

Error structure * autocorrelative errors

library(nlme)

# # With RE (probably can be done)
# # incompatible formulas for groups in 'random' and 'correlation'
# temp_model_ar <- wigley_bmspectra_model_df %>% 
#   mutate(group = str_c(survey_area, "_", season)) %>% 
#   nlme::lme(
#   fixed = b ~ survey_area * season * scale(bot_temp_anom), 
#   random = ~ 1 | est_year, 
#   correlation = corAR1(form = ~ est_year | survey_area/season), 
#   data = .
# ) 


# Just drop the random effect part for now
temp_model_ar <- nlme::gls(
  b ~ survey_area * season * scale(bot_temp_anom), 
  data = wigley_bmspectra_model_df, 
  correlation = corAR1(form = ~ est_year | survey_area/season)
)


# # confidence intervals on phi
# intervals(temp_model_ar)$corStruct


# # check the model
# check_model(temp_model_ar)
# plot(check_collinearity(temp_model_ar))


# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    temp_model_ar, 
    terms = ~ bot_temp_anom * survey_area * season))   %>% 
  mutate(
    survey_area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = temp_model_ar,
  ~survey_area * season,
  var = "bot_temp_anom")


# Drop effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = temp_model_ar,
  specs =  ~ survey_area * season,
  var = "bot_temp_anom") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, survey_area) %>%
  summarise(
    min_temp = min(bot_temp_anom)-2,
    max_temp = max(bot_temp_anom)+2,
    .groups = "drop")



# Plot over observed data
local_btemp_anoms <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(bot_temp_anom, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(survey_area~., scales = "free") +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Exponent of ISD",
       title = "Local Bottom Temperature Anomalies & area * season Relationships",
       x = "Local Seasonal Bottom Temperature Anomaly")
local_btemp_anoms

Based on this model I would say there is evidence of a warming effect in: GOM-Spring, GB-Spring, GB-Fall, & SNE-Fall.

Total Landings

The landings information is the newer GARFO data obtained from Andy Beet. It was delivered as containing only finfish, I have removed freshwater species and blue crabs. It should largely reflect marine finfish species at this point.

# From Andy Beet - Crabs removed
landings_annual <- read_csv(here::here("Data/processed/BEET_GARFO_regional_finfish_landings.csv")) %>% 
  select(survey_area,
         est_year = year,
         total_weight_lb = total_live_lb)

shelf_landings_annual <-  mutate(landings_annual, survey_area = "Northeast Shelf") %>% 
    group_by(est_year, survey_area) %>% 
    summarise(total_weight_lb = sum(total_weight_lb, na.rm = T))

# Plot annual totals
bind_rows(list(
  shelf_landings_annual,
  landings_annual)) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels_long_all)) %>% 
  filter(est_year <= 2019) %>% 
  ggplot(aes(est_year, total_weight_lb/1e6)) +
  geom_col(alpha = 0.5) +
  geom_ma(aes(color = "5-Year Smooth"), size = 1, n = 5, ma_fun = SMA, linetype = 1) +
  scale_color_manual(values = "orange") +
  facet_wrap(~survey_area, ncol = 1, scales = "free") +
  scale_y_continuous(labels = label_number(suffix = "M")) +
  labs(y = "Total Landings (live lb.)", x = "Year", color = "Moving Average")

Landings Autoregressive model

I wanted to follow on the bottom temperature anomalies model here for landings:

library(nlme)

# # With RE (probably can be done)
# # incompatible formulas for groups in 'random' and 'correlation'
# temp_model_ar <- wigley_bmspectra_model_df %>% 
#   mutate(group = str_c(survey_area, "_", season)) %>% 
#   nlme::lme(
#   fixed = b ~ survey_area * season * scale(bot_temp_anom), 
#   random = ~ 1 | est_year, 
#   correlation = corAR1(form = ~ est_year | survey_area/season), 
#   data = .
# ) 


# Just drop the random effect part for now
f_model_ar <- nlme::gls(
  b ~ survey_area * season + survey_area * scale(log(total_weight_lb)), 
  data = drop_na(wigley_bmspectra_model_df, total_weight_lb),
  correlation = corAR1(form = ~ est_year | survey_area/season)
)


# # confidence intervals on phi
# intervals(temp_model_ar)$corStruct


# # check the model
# check_model(f_model_ar)
# plot(check_collinearity(f_model_ar))


# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    f_model_ar, 
    terms = ~ total_weight_lb * survey_area * season))   %>% 
  mutate(
    survey_area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = f_model_ar,
  ~survey_area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite")


# Drop effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = f_model_ar,
  specs =  ~ survey_area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, survey_area) %>%
  summarise(
    min_f = min(total_weight_lb)-2,
    max_f = max(total_weight_lb)+2,
    .groups = "drop")



# Plot over observed data
landings_ar_plot <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_f) == F,
         (x > max_f) == F) %>%
  left_join(select(mod2_emtrend, survey_area,non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(survey_area~., scales = "free") +
  scale_color_gmri() +
  scale_x_log10(labels = label_log(base = 10), limits = 10^c(0,10)) +
  labs(y = "Exponent of ISD",
       title = "Total Landings (lb.)",
       x = "Local Seasonal Bottom Temperature Anomaly")
landings_ar_plot

# temp model is better
# AIC(temp_model_ar)
# AIC(f_model_ar)

Zooplankton Community Indices

The ecodata package has changed the zooplankton indices around that it carries. It now contains “regime” metrics for 26 of the main zooplankton taxa.

https://noaa-edab.github.io/tech-doc/zoo_abundance_anom.html

zoo_regimes <- ecodata::zoo_regime
# zoo_anoms <- ecodata::zoo_abundance_anom
zoo_regimes %>% 
  filter(str_detect(Var, "calfin|chaeto|pseudo|oith")) %>% 
  mutate(EPU = factor(EPU, levels = c("GOM", "GB", "MAB"))) %>% 
ggplot(aes(Time, Value)) +
  geom_point(
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(linetype = "5-Year Smooth"),
    size = 1, n = 5, ma_fun = SMA) +
  #scale_color_gmri() +
  facet_grid(EPU~Var) +
  labs()

---

Driver Correlations

Modeling Driver Relationships

This section is likely to change

The following model has been used to determine the relationship between landings and seasonal bottom temperature on b.

The Northeast Shelf as a whole is not included as its own region in these models, but is plotted below with points just for visual context.

lmer(
  b ~ survey_area*season*scale(roll_temp) + survey_area * scale(log(total_weight_lb)) + (1|est_year), 
  data = wtb_model_df)

#### Model 2: Temperature & Fishing ####


# # Drop NA's
# wtb_model_df <- 
#   drop_na(wigley_bmspectra_model_df, total_weight_lb, bot_temp, b) %>%
#   rename(area = survey_area) %>% 
#   left_join(area_df) %>% 
  
  
wtb_model_df <- drop_na(wigley_bmspectra_model_df, bot_temp, b) %>% 
  # Do rolling BT within a season * region
  group_by(survey_area, season) %>%
  mutate(
    roll_temp = zoo::rollapply(bot_temp, 5, mean, na.rm = T, align = "right",  fill = NA)) %>% 
  ungroup() %>% 
  mutate(
    yr_num = as.numeric(est_year),
    yr_fac = factor(est_year),
    survey_area = factor(survey_area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall")),
    landings = total_weight_lb,
    yr_seas = str_c(season, est_year))



# Make the model
mod2 <- lmer(
  b ~ survey_area*season*scale(roll_temp) + survey_area*scale(log(total_weight_lb)) + (1|est_year), 
  data = wtb_model_df)


# summary(mod2)
# check_model(mod2)

Temperature Relationship

If we follow yearly region * season combinations as our units of comparison the relationship between seasonal bottom temperature and the size distribution falls apart.

# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    mod2, 
    terms = ~ roll_temp*survey_area*season))   %>% 
  mutate(
    survey_area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = mod2,
  ~survey_area * season,
  var = "roll_temp")


# Drop effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = mod2,
  specs =  ~ survey_area*season,
  var = "roll_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))



# Limit temperature plotting range
actual_range <- wtb_model_df %>% 
  group_by(season, survey_area) %>% 
  summarise(min_temp = min(bot_temp)-2,
            max_temp = max(bot_temp)+2,
            .groups = "drop")



# Plot over observed data
btemp_local <- mod2_preds %>% 
  left_join(actual_range) %>% 
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>% 
  left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = survey_area), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wtb_model_df,
    aes(roll_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(survey_area~., scales = "free") +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Exponent of ISD",
       title = "Within area*season Relationship",
       x = "5-Year Rolling Average Bottom Temperature")
btemp_local

A spring spectra is not any more/less steep if spring is hotter or colder than usual. However, Fall is consistently steeper than spring, and fall is always hotter than spring.

Landings Relationship

Since landings are annual there won’t be any seasonal interaction, only an intercept adjustment for the survey season.

# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    mod2, 
    terms = ~ total_weight_lb * survey_area * season))   %>% 
  mutate(
    survey_area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = mod2,
  ~survey_area * season,
  var = "total_weight_lb")
#summary(trend_df, infer= c(T,T))


# Drop effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = mod2,
  specs =  ~ survey_area*season,
  var = "log(total_weight_lb)") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T),
    group = str_c(survey_area, "X", season))




# Plot over observed data
mod2_preds %>% 
  left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wtb_model_df,
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(survey_area~., scales = "free") +
  scale_color_gmri() +
  scale_x_continuous(transform = "log10", labels = label_log(base = 10)) +
  labs(y = "Body Mass Spectra Slope (b)",
       title = "Wigley Species, Body Mass Spectra",
       x = "Total Landings (lb.)")

Landings + Temperature Modeling Considerations

With the old landings data we had trouble determining attribution from the models with both temperature and landings simultaneously due to strong collinearities. This may no longer be the case, lets look.

The bigger issue looks to be the relationship between survey season + region and bottom temperature. There are other phenological dynamics (spring plankton bloom, seasonal migrations, spawning, stratification etc.) that happen seasonally that I am hesitant to use temperature alone without controlling for season.

If we check variance inflation factors for the model we currently have nearly all predictors have very high VIF values, beyond what is the typical cutoff.

Which is unfortunate, because it seems like its a pretty reasonable model otherwise.

# Drop effect fits that are non-significant  ###


# summary(mod2)
# summary(temp_landings_mod)
check_model(mod2)

Driver Model Summary

# mod2 %>% 
#   gtsummary::tbl_regression()   %>% 
#   add_global_p(keep = T) %>% 
#   bold_p(t = 0.05)  %>%
#   bold_labels() %>%
#   italicize_levels() %>% 
#   as_gt() %>% 
#   tab_header(title = "Wigley Community - Bodymass Distribution & Drivers")

Side Stories

Everything below is digging into different nuances of the approach or the community composition. These things are of secondary importance and will be moved out.

Spiny Dogfish Migrations

Spiny dogfish constitute a majority fraction of biomass in the MAB during the Spring survey, but are subsequently absent during the fall season.

Their numbers have also been increasing over time and they are unquestionably the dominant “large” species caught in the survey.

Their growing abundance, particularly their spring abundance in the Mid-Atlantic Bight, acts to shallow the spectra as they are one of the larger species routinely sampled by the trawl survey.

If we remove them from the estimation and re-run the slopes this is what changes.

Once removed we see that the Mid-Atlantic Bight is no longer becoming less steep over time, with no long-term trend. Spectra in Georges Bank without the inclusion of spiny dogfish are now declining in both seasons.

# # Run MAB with and without spiny dogfish
# # Run the year, season, area fits for the filtered data
# no_dogs_bmspectra <- group_binspecies_spectra(
#   ss_input = filter(wigley_in, comname != "spiny dogfish"),
#   grouping_vars = c("est_year", "season", "survey_area"),
#   abundance_vals = "numlen_adj",
#   length_vals = "length_cm",
#   use_weight = TRUE,
#   isd_xmin = 1,
#   global_min = TRUE,
#   isd_xmax = NULL,
#   global_max = FALSE,
#   bin_width = 1,
#   vdiff = 2)

# # Save those
# write_csv(
#   no_dogs_bmspectra,
#   here::here("Data/processed/wigley_species_nodogfish_bmspectra.csv"))

# Read them in, do some plotting
no_dogs_bmspectra <- read_csv(
  here::here("Data/processed/wigley_species_nodogfish_bmspectra.csv")) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels))


# Format a little
no_dogs_bmspectra <- no_dogs_bmspectra %>% 
  mutate(est_year = as.numeric(est_year))


# Join them together
dfish_results <- bind_rows(
   list(
     "Wigley Full" = wigley_bmspectra,
     "Wigley w/o Spiny Dogfish" =  no_dogs_bmspectra), 
   .id = "community") %>% 
  mutate(metric = "bodymass_spectra") %>% 
  mutate(
    survey_area = fct_relevel(survey_area, area_levels))




# Get trends for no dogfish
nodfish_mod <- lmer(
  formula = b ~ est_year * survey_area * season + (1|est_year),
  data = no_dogs_bmspectra)
bmspectra_mod_wig <- lmer(
  formula = b ~ est_year * survey_area * season + (1|est_year),
  data = wigley_bmspectra)

# Put predictions
dogfish_sens_predictions <- bind_rows(list(
  "bodymass_spectra-Wigley Full" = get_preds_and_trendsignif(bmspectra_mod_wig),
  "bodymass_spectra-Wigley w/o Spiny Dogfish" = get_preds_and_trendsignif(nodfish_mod)), 
  .id = "model_id") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric))



# Plot the differences
ggplot() +
  geom_ribbon(
    data = filter(dogfish_sens_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    data = dfish_results,
    aes(est_year, b, color = season, 
        alpha = I(if_else(survey_area %in% c("MAB", "GB"), 1, 0.2))),
    size = 0.8) +
  geom_line(
    data = filter(dogfish_sens_predictions, non_zero == TRUE),
    aes(x, predicted, color = season,
        alpha = I(if_else(survey_area %in% c("MAB", "GB"), 1, 0.2))),
    linewidth = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(survey_area ~ metric*community, scales = "free") +
  labs(
    x = "Year",
    y = "Exponent of Size Spectra",
    title = "Shift in Spectra Trend(s) with Exclusion of Spiny Dogfish",
    subtitle = "Spring = Dogfish in MAB, Fall = Dogfish in GB")

Spiny Dogfish Seasonal Dstributions

There is literature documenting a seasonal shift in center of biomass Northward in autumn, and Southward in the spring:

Sagarese et a. 2015.

And another citation documenting their propensity to migrate and range large distances Carlson et al. 20

And some documentation of sexual segregation and habitat usage: Janne B. Haugen et al. 2017

Story Thoughts

Below is a table of some papers operating in the same area or around the same topic that are particularly relevant:

shelf_papers <- tribble(
  ~"Author", ~"Year", ~"Conjecture",
  "Friedland et al.", 2024, "Body size has declined despite more abundance & biomass. Competition is preventing larger sizes.",
  "Krumsick", 2020, "Empirical slopes from trawl survey in labrador sea were steeper than theoretical from nitrogen stable isotope ratios"
)


shelf_papers %>% gt()

Author	Year	Conjecture
Friedland et al.	2024	Body size has declined despite more abundance & biomass. Competition is preventing larger sizes.
Krumsick	2020	Empirical slopes from trawl survey in labrador sea were steeper than theoretical from nitrogen stable isotope ratios

There are a couple core ideas about what changes have occurred or have been ocurring over the Northeast Shelf with relevance to the community that is sampled by the trawl survey: